Chapter 6 Beta diversity

load("data/data.Rdata")

beta_q0n <- genome_counts_filt %>%
  column_to_rownames(., "genome") %>%
  filter(rowSums(. != 0, na.rm = TRUE) > 0) %>%
  dplyr::select_if(~!all(. == 0)) %>%
  hillpair(., q = 0)

beta_q1n <- genome_counts_filt %>%
  column_to_rownames(., "genome") %>%
  filter(rowSums(. != 0, na.rm = TRUE) > 0) %>%
  dplyr::select_if(~!all(. == 0)) %>%
  hillpair(., q = 1)

beta_q1p <- genome_counts_filt %>%
  column_to_rownames(., "genome") %>%
  filter(rowSums(. != 0, na.rm = TRUE) > 0) %>%
  dplyr::select_if(~!all(. == 0)) %>%
  hillpair(., q = 1, tree = genome_tree)

beta_q1f <- genome_counts_filt %>%
  column_to_rownames(., "genome") %>%
  filter(rowSums(. != 0, na.rm = TRUE) > 0) %>%
  dplyr::select_if(~!all(. == 0)) %>%
  hillpair(., q = 1, dist = dist)

6.1 location

6.1.1 Richness diversity

betadisper(beta_q0n$S, sample_metadata$location) %>% permutest(., pairwise=TRUE) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)    
Groups     5 0.51463 0.102926 11.076    999  0.001 ***
Residuals 86 0.79917 0.009293                         
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
               Aruba     Brazil  CaboVerde      Spain    Denmark Malaysia
Aruba                1.0000e-03 2.0000e-03 1.0000e-03 1.0000e-03    0.010
Brazil    8.1559e-05            1.1500e-01 7.1000e-02 9.1900e-01    0.135
CaboVerde 6.2922e-04 1.1254e-01            1.0000e-03 2.6000e-02    0.930
Spain     1.7363e-12 7.3505e-02 3.5940e-05            2.7000e-02    0.002
Denmark   1.3206e-07 9.0556e-01 2.9720e-02 2.3526e-02               0.059
Malaysia  9.6128e-03 1.3730e-01 9.2726e-01 4.6570e-04 5.9273e-02         

adonis2(beta_q0n$S ~ location, 
        data = sample_metadata %>% arrange(match(sample,labels(beta_q0n$S))), 
        permutations = 999) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
location	5	6.57607	0.2442751	5.559605	0.001
Residual	86	20.34468	0.7557249	NA	NA
Total	91	26.92075	1.0000000	NA	NA

pairwise.adonis(beta_q0n$S, sample_metadata$location, perm = 999)

                   pairs Df SumsOfSqs  F.Model         R2 p.value p.adjusted sig
1        Aruba vs Brazil  1 2.0520609 7.280171 0.19528265   0.001      0.015   .
2     Aruba vs CaboVerde  1 1.0026264 3.213494 0.09975614   0.001      0.015   .
3       Aruba vs Denmark  1 1.9513181 7.018855 0.19486614   0.001      0.015   .
4      Aruba vs Malaysia  1 1.1650328 3.678697 0.11257172   0.001      0.015   .
5         Aruba vs Spain  1 2.0433431 8.167535 0.21974917   0.001      0.015   .
6    Brazil vs CaboVerde  1 2.1090004 9.164611 0.24013374   0.001      0.015   .
7      Brazil vs Denmark  1 0.5703370 2.907951 0.09113563   0.001      0.015   .
8     Brazil vs Malaysia  1 0.9283278 3.953419 0.11996993   0.001      0.015   .
9        Brazil vs Spain  1 0.6343957 3.769479 0.11503016   0.001      0.015   .
10  CaboVerde vs Denmark  1 1.8099512 8.070070 0.22373314   0.001      0.015   .
11 CaboVerde vs Malaysia  1 1.2143821 4.593887 0.14094322   0.001      0.015   .
12    CaboVerde vs Spain  1 1.8141430 9.281722 0.24896172   0.001      0.015   .
13   Denmark vs Malaysia  1 1.0124763 4.418610 0.13629855   0.001      0.015   .
14      Denmark vs Spain  1 0.5305346 3.310771 0.10573905   0.001      0.015   .
15     Malaysia vs Spain  1 0.8450695 4.218747 0.13094076   0.001      0.015   .

#pdf("figures/beta_q0_loca.pdf",width=9, height=5)
beta_q0n$S %>%
  vegan::metaMDS(., trymax = 500, k = 2, trace=0) %>%
  vegan::scores() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(sample_metadata, by = "sample") %>%
  group_by(location) %>%
  mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
  mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
  ungroup() %>%
  ggplot(aes(x = NMDS1, y = NMDS2, color = location, fill = location)) +
    geom_point(size = 4) +
    #   stat_ellipse(aes(color = beta_q1n_nmds$Groups))+
    geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2), alpha = 0.9, show.legend = FALSE) +
    scale_color_manual(values = location_colors)+
    theme_classic() +
    theme(
      axis.text.x = element_text(size = 12),
      axis.text.y = element_text(size = 12),
      axis.title = element_text(size = 20, face = "bold"),
      axis.text = element_text(face = "bold", size = 18),
      panel.background = element_blank(),
      axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
      legend.text = element_text(size = 16),
      legend.title = element_text(size = 18),
      legend.position = "right", legend.box = "vertical"
    ) 

#dev.off()

6.1.2 Neutral diversity

betadisper(beta_q1n$S, sample_metadata$location) %>% permutest(., pairwise=TRUE) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df Sum Sq  Mean Sq      F N.Perm Pr(>F)   
Groups     5 0.2728 0.054560 3.7079    999  0.003 **
Residuals 86 1.2654 0.014715                        
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
               Aruba     Brazil  CaboVerde      Spain    Denmark Malaysia
Aruba                9.0000e-03 1.0000e-03 1.0000e-03 2.0000e-03    0.002
Brazil    9.2383e-03            7.7300e-01 2.0900e-01 6.5800e-01    0.872
CaboVerde 7.7521e-04 7.7033e-01            2.8900e-01 8.7700e-01    0.905
Spain     1.6646e-05 2.2327e-01 3.0424e-01            3.7300e-01    0.262
Denmark   3.3254e-04 6.6322e-01 8.7642e-01 3.7291e-01               0.795
Malaysia  1.1422e-03 8.5989e-01 8.9763e-01 2.4754e-01 7.7435e-01         

adonis2(beta_q1n$S ~ location, 
        data = sample_metadata %>% arrange(match(sample,labels(beta_q1n$S))), 
        permutations = 999) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
location	5	5.709562	0.221304	4.88821	0.001
Residual	86	20.090068	0.778696	NA	NA
Total	91	25.799630	1.000000	NA	NA

pairwise.adonis(beta_q1n$S, sample_metadata$location, perm = 999)

                   pairs Df SumsOfSqs  F.Model         R2 p.value p.adjusted sig
1        Aruba vs Brazil  1 1.3739907 4.789154 0.13766227   0.001      0.015   .
2     Aruba vs CaboVerde  1 1.0636589 3.806728 0.11603497   0.001      0.015   .
3       Aruba vs Denmark  1 1.6037782 5.813780 0.16699651   0.001      0.015   .
4      Aruba vs Malaysia  1 1.2584980 4.466640 0.13346545   0.001      0.015   .
5         Aruba vs Spain  1 1.4382789 5.558030 0.16083181   0.001      0.015   .
6    Brazil vs CaboVerde  1 1.5901673 7.027154 0.19505161   0.001      0.015   .
7      Brazil vs Denmark  1 0.6954159 3.122217 0.09719804   0.004      0.060    
8     Brazil vs Malaysia  1 0.7313843 3.199010 0.09935119   0.001      0.015   .
9        Brazil vs Spain  1 0.4007935 1.948927 0.06297236   0.033      0.495    
10  CaboVerde vs Denmark  1 1.8198218 8.556140 0.23405480   0.001      0.015   .
11 CaboVerde vs Malaysia  1 1.3171418 6.019861 0.17695136   0.001      0.015   .
12    CaboVerde vs Spain  1 1.5415454 7.905417 0.22017337   0.001      0.015   .
13   Denmark vs Malaysia  1 1.0691578 4.970172 0.15074753   0.001      0.015   .
14      Denmark vs Spain  1 0.5766622 3.014217 0.09718823   0.002      0.030   .
15     Malaysia vs Spain  1 0.6446972 3.265586 0.10444666   0.002      0.030   .

#pdf("figures/beta_q1n_loca.pdf",width=9, height=5)
beta_q1n$S %>%
  vegan::metaMDS(., trymax = 500, k = 2, trace=0) %>%
  vegan::scores() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(sample_metadata, by = "sample") %>%
  group_by(location) %>%
  mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
  mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
  ungroup() %>%
  ggplot(aes(x = NMDS1, y = NMDS2, color = location, fill = location)) +
    geom_point(size = 4) +
    #   stat_ellipse(aes(color = beta_q1n_nmds$Groups))+
    geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2), alpha = 0.9, show.legend = FALSE) +
    scale_color_manual(values = location_colors)+
    theme_classic() +
    theme(
      axis.text.x = element_text(size = 12),
      axis.text.y = element_text(size = 12),
      axis.title = element_text(size = 20, face = "bold"),
      axis.text = element_text(face = "bold", size = 18),
      panel.background = element_blank(),
      axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
      legend.text = element_text(size = 16),
      legend.title = element_text(size = 18),
      legend.position = "right", legend.box = "vertical"
    ) 

#dev.off()

6.1.3 Phylogenetic diversity

betadisper(beta_q1p$S, sample_metadata$location) %>% permutest(., pairwise=TRUE) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df Sum Sq Mean Sq      F N.Perm Pr(>F)  
Groups     5 0.2152 0.04304 2.1047    999  0.081 .
Residuals 86 1.7587 0.02045                       
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
               Aruba     Brazil  CaboVerde      Spain    Denmark Malaysia
Aruba                0.06000000 0.04600000 0.00200000 0.11100000    0.005
Brazil    0.06236518            0.87300000 0.32400000 0.65100000    0.712
CaboVerde 0.05279004 0.86912737            0.24300000 0.76500000    0.545
Spain     0.00060645 0.32391452 0.20693940            0.12100000    0.514
Denmark   0.09811567 0.63811338 0.73976492 0.10108967               0.324
Malaysia  0.00519311 0.67845753 0.52055144 0.50696000 0.31212601         

adonis2(beta_q1p$S ~ location, 
        data = sample_metadata %>% arrange(match(sample,labels(beta_q1p$S))), 
        permutations = 999) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
location	5	2.690791	0.2683828	6.309563	0.001
Residual	86	7.335152	0.7316172	NA	NA
Total	91	10.025944	1.0000000	NA	NA

pairwise.adonis(beta_q1p$S, sample_metadata$location, perm = 999)

                   pairs Df  SumsOfSqs    F.Model         R2 p.value p.adjusted sig
1        Aruba vs Brazil  1 0.64714035  6.1133065 0.16928127   0.001      0.015   .
2     Aruba vs CaboVerde  1 0.56229778  5.3504131 0.15575979   0.001      0.015   .
3       Aruba vs Denmark  1 1.03758657  9.5591181 0.24790811   0.001      0.015   .
4      Aruba vs Malaysia  1 0.48047063  5.0352414 0.14794199   0.001      0.015   .
5         Aruba vs Spain  1 0.60047832  6.7919744 0.18976250   0.001      0.015   .
6    Brazil vs CaboVerde  1 0.70961443  8.1023431 0.21837821   0.001      0.015   .
7      Brazil vs Denmark  1 0.32742702  3.5968626 0.11034383   0.025      0.375    
8     Brazil vs Malaysia  1 0.19570506  2.5119802 0.07971509   0.028      0.420    
9        Brazil vs Spain  1 0.06456161  0.9106384 0.03044530   0.456      1.000    
10  CaboVerde vs Denmark  1 1.07120011 11.9405505 0.29895808   0.001      0.015   .
11 CaboVerde vs Malaysia  1 0.50032095  6.5728056 0.19011490   0.002      0.030   .
12    CaboVerde vs Spain  1 0.82104026 11.9236999 0.29866220   0.001      0.015   .
13   Denmark vs Malaysia  1 0.54258755  6.8084757 0.19559821   0.001      0.015   .
14      Denmark vs Spain  1 0.24626102  3.3999417 0.10827860   0.019      0.285    
15     Malaysia vs Spain  1 0.26932060  4.5771932 0.14050299   0.001      0.015   .

#pdf("figures/beta_q1p_loca.pdf",width=9, height=5)
beta_q1p$S %>%
  vegan::metaMDS(., trymax = 500, k = 2, trace=0) %>%
  vegan::scores() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(sample_metadata, by = "sample") %>%
  group_by(location) %>%
  mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
  mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
  ungroup() %>%
  ggplot(aes(x = NMDS1, y = NMDS2, color = location, fill = location)) +
    geom_point(size = 4) +
    #   stat_ellipse(aes(color = beta_q1n_nmds$Groups))+
    geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2), alpha = 0.9, show.legend = FALSE) +
    scale_color_manual(values = location_colors)+
    theme_classic() +
    theme(
      axis.text.x = element_text(size = 12),
      axis.text.y = element_text(size = 12),
      axis.title = element_text(size = 20, face = "bold"),
      axis.text = element_text(face = "bold", size = 18),
      panel.background = element_blank(),
      axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
      legend.text = element_text(size = 16),
      legend.title = element_text(size = 18),
      legend.position = "right", legend.box = "vertical"
    ) 

#dev.off()

6.1.4 Functional diversity

betadisper(beta_q1f$S, sample_metadata$location) %>% permutest(., pairwise=TRUE) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df Sum Sq  Mean Sq      F N.Perm Pr(>F)    
Groups     5 1.5016 0.300325 9.1732    999  0.001 ***
Residuals 86 2.8156 0.032739                         
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
               Aruba     Brazil  CaboVerde      Spain    Denmark Malaysia
Aruba                2.0000e-03 5.9800e-01 1.0000e-03 1.0000e-03    0.475
Brazil    8.6955e-04            4.0000e-03 6.9600e-01 2.4600e-01    0.002
CaboVerde 5.8936e-01 9.7727e-04            3.0000e-03 1.0000e-03    0.893
Spain     1.5125e-03 6.8358e-01 1.5804e-03            3.7000e-02    0.002
Denmark   1.7539e-04 2.3957e-01 9.2505e-05 3.3658e-02               0.001
Malaysia  4.9615e-01 1.0071e-03 8.8292e-01 1.5959e-03 7.4748e-05         

adonis2(beta_q1f$S ~ location, 
        data = sample_metadata %>% arrange(match(sample,labels(beta_q1f$S))), 
        permutations = 999) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
location	5	8.815599	0.5844961	24.19552	0.001
Residual	86	6.266793	0.4155039	NA	NA
Total	91	15.082391	1.0000000	NA	NA

pairwise.adonis(beta_q1f$S, sample_metadata$location, perm = 999)

                   pairs Df   SumsOfSqs    F.Model           R2 p.value p.adjusted sig
1        Aruba vs Brazil  1  3.77109355 40.6261637  0.575228238   0.001      0.015   .
2     Aruba vs CaboVerde  1  0.81599223  5.5801053  0.161367505   0.008      0.120    
3       Aruba vs Denmark  1  4.54372380 50.8948430  0.637022880   0.001      0.015   .
4      Aruba vs Malaysia  1  0.87672195  6.0397240  0.172367910   0.007      0.105    
5         Aruba vs Spain  1  3.49955224 36.9984243  0.560595570   0.001      0.015   .
6    Brazil vs CaboVerde  1  1.79573937 27.7696018  0.489163230   0.001      0.015   .
7      Brazil vs Denmark  1  0.27026868 35.0545409  0.547260825   0.001      0.015   .
8     Brazil vs Malaysia  1  1.69391038 26.6369446  0.478763612   0.001      0.015   .
9        Brazil vs Spain  1 -0.01258568 -0.9666596 -0.034482498   0.946      1.000    
10  CaboVerde vs Denmark  1  2.85656348 47.6207158  0.629731090   0.001      0.015   .
11 CaboVerde vs Malaysia  1 -0.02164970 -0.1836835 -0.006603446   0.984      1.000    
12    CaboVerde vs Spain  1  1.56212224 23.8545897  0.460028512   0.001      0.015   .
13   Denmark vs Malaysia  1  2.82610262 48.0023359  0.631590271   0.001      0.015   .
14      Denmark vs Spain  1  0.38029019 58.5480612  0.676480332   0.001      0.015   .
15     Malaysia vs Spain  1  1.47580731 22.9256154  0.450178466   0.001      0.015   .

#pdf("figures/beta_q1f_loca.pdf",width=9, height=5)
beta_q1f$S %>%
  vegan::metaMDS(., trymax = 500, k = 2, trace=0) %>%
  vegan::scores() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(sample_metadata, by = "sample") %>%
  group_by(location) %>%
  mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
  mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
  ungroup() %>%
  ggplot(aes(x = NMDS1, y = NMDS2, color = location, fill = location)) +
    geom_point(size = 4) +
    #   stat_ellipse(aes(color = beta_q1n_nmds$Groups))+
    geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2), alpha = 0.9, show.legend = FALSE) +
    scale_color_manual(values = location_colors)+
    theme_classic() +
    theme(
      axis.text.x = element_text(size = 12),
      axis.text.y = element_text(size = 12),
      axis.title = element_text(size = 20, face = "bold"),
      axis.text = element_text(face = "bold", size = 18),
      panel.background = element_blank(),
      axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
      legend.text = element_text(size = 16),
      legend.title = element_text(size = 18),
      legend.position = "right", legend.box = "vertical"
    ) 

#dev.off()

6.2 Behaviour

6.2.1 Richness diversity

betadisper(beta_q0n$S, sample_metadata$origin) %>% permutest(., pairwise=TRUE) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     1 0.00411 0.0041147 0.2442    999  0.602
Residuals 90 1.51634 0.0168483                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
         Domestic Feral
Domestic            0.6
Feral     0.62238      

adonis2(beta_q0n$S ~ origin, 
        data = sample_metadata %>% arrange(match(sample,labels(beta_q0n$S))), 
        permutations = 999,
        strata = sample_metadata %>% arrange(match(sample,labels(beta_q1n$S))) %>% pull(location)) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
origin	1	0.4407951	0.0163738	1.498173	0.194
Residual	90	26.4799587	0.9836262	NA	NA
Total	91	26.9207538	1.0000000	NA	NA

pairwise.adonis(beta_q0n$S, sample_metadata$origin, perm = 999)

              pairs Df SumsOfSqs  F.Model        R2 p.value p.adjusted sig
1 Domestic vs Feral  1 0.4407951 1.498173 0.0163738   0.084      0.084    

beta_q0n$S %>%
  vegan::metaMDS(., trymax = 500, k = 2, trace=0) %>%
  vegan::scores() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(sample_metadata, by = join_by(sample == sample)) %>%
  group_by(origin) %>%
  mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
  mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
  ungroup() %>%
  ggplot(aes(x = NMDS1, y = NMDS2, color = origin, fill = origin)) +
    geom_point(size = 4) +
    scale_color_manual(values = origin_colors) +
    geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2), alpha = 0.9, show.legend = FALSE) +
    theme_classic() +
    theme(
      axis.text.x = element_text(size = 12),
      axis.text.y = element_text(size = 12),
      axis.title = element_text(size = 20, face = "bold"),
      axis.text = element_text(face = "bold", size = 18),
      panel.background = element_blank(),
      axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
      legend.text = element_text(size = 16),
      legend.title = element_text(size = 18),
      legend.position = "right", legend.box = "vertical"
    ) 

6.2.2 Neutral diversity

betadisper(beta_q1n$S, sample_metadata$origin) %>% permutest(., pairwise=TRUE) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     1 0.0080 0.0080048 0.5581    999  0.464
Residuals 90 1.2909 0.0143435                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
         Domestic Feral
Domestic          0.465
Feral     0.45698      

adonis2(beta_q1n$S ~ origin, 
        data = sample_metadata %>% arrange(match(sample,labels(beta_q1n$S))), 
        permutations = 999,
        strata = sample_metadata %>% arrange(match(sample,labels(beta_q1n$S))) %>% pull(location)) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
origin	1	0.386668	0.01498735	1.369385	0.302
Residual	90	25.412962	0.98501265	NA	NA
Total	91	25.799630	1.00000000	NA	NA

pairwise.adonis(beta_q1n$S, sample_metadata$origin, perm = 999)

              pairs Df SumsOfSqs  F.Model         R2 p.value p.adjusted sig
1 Domestic vs Feral  1  0.386668 1.369385 0.01498735   0.142      0.142    

pdf("figures/beta_q1n_behaviour.pdf",width=9, height=5)
beta_q1n$S %>%
  vegan::metaMDS(., trymax = 500, k = 2, trace=0) %>%
  vegan::scores() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(sample_metadata, by = join_by(sample == sample)) %>%
  group_by(origin) %>%
  mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
  mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
  ungroup() %>%
  ggplot(aes(x = NMDS1, y = NMDS2, color = origin, fill = origin)) +
    geom_point(size = 4) +
    scale_color_manual(values = origin_colors) +
    geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2), alpha = 0.9, show.legend = FALSE) +
    theme_classic() +
    theme(
      axis.text.x = element_text(size = 12),
      axis.text.y = element_text(size = 12),
      axis.title = element_text(size = 20, face = "bold"),
      axis.text = element_text(face = "bold", size = 18),
      panel.background = element_blank(),
      axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
      legend.text = element_text(size = 16),
      legend.title = element_text(size = 18),
      legend.position = "right", legend.box = "vertical"
    ) 
dev.off()

quartz_off_screen 
                2 

6.2.3 Phylogenetic diversity

betadisper(beta_q1p$S, sample_metadata$origin) %>% permutest(., pairwise=TRUE) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)
Groups     1 0.01229 0.012295 0.5442    999  0.451
Residuals 90 2.03337 0.022593                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
         Domestic Feral
Domestic          0.452
Feral     0.46262      

adonis2(beta_q1p$S ~ origin, 
        data = sample_metadata %>% arrange(match(sample,labels(beta_q1p$S))), 
        permutations = 999,
        strata = sample_metadata %>% arrange(match(sample,labels(beta_q1n$S))) %>% pull(location)) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
origin	1	0.2191073	0.02185403	2.010807	0.183
Residual	90	9.8068363	0.97814597	NA	NA
Total	91	10.0259435	1.00000000	NA	NA

pairwise.adonis(beta_q1p$S, sample_metadata$origin, perm = 999)

              pairs Df SumsOfSqs  F.Model         R2 p.value p.adjusted sig
1 Domestic vs Feral  1 0.2191073 2.010807 0.02185403   0.073      0.073    

beta_q1p$S %>%
  vegan::metaMDS(., trymax = 500, k = 2, trace=0) %>%
  vegan::scores() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(sample_metadata, by = join_by(sample == sample)) %>%
  group_by(origin) %>%
  mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
  mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
  ungroup() %>%
  ggplot(aes(x = NMDS1, y = NMDS2, color = origin, fill = origin)) +
    geom_point(size = 4) +
    scale_color_manual(values = origin_colors) +
    geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2), alpha = 0.9, show.legend = FALSE) +
    theme_classic() +
    theme(
      axis.text.x = element_text(size = 12),
      axis.text.y = element_text(size = 12),
      axis.title = element_text(size = 20, face = "bold"),
      axis.text = element_text(face = "bold", size = 18),
      panel.background = element_blank(),
      axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
      legend.text = element_text(size = 16),
      legend.title = element_text(size = 18),
      legend.position = "right", legend.box = "vertical"
    ) 

6.2.4 Functional diversity

betadisper(beta_q1f$S, sample_metadata$origin) %>% permutest(., pairwise=TRUE) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df Sum Sq  Mean Sq      F N.Perm Pr(>F)  
Groups     1 0.2486 0.248604 3.3502    999  0.058 .
Residuals 90 6.6786 0.074206                       
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
         Domestic Feral
Domestic          0.061
Feral    0.070508      

adonis2(beta_q1f$S ~ origin, 
        data = sample_metadata %>% arrange(match(sample,labels(beta_q1f$S))), 
        permutations = 999,
        strata = sample_metadata %>% arrange(match(sample,labels(beta_q1n$S))) %>% pull(location)) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
origin	1	0.3816603	0.02530502	2.336579	0.198
Residual	90	14.7007312	0.97469498	NA	NA
Total	91	15.0823914	1.00000000	NA	NA

pairwise.adonis(beta_q1f$S, sample_metadata$origin, perm = 999)

              pairs Df SumsOfSqs  F.Model         R2 p.value p.adjusted sig
1 Domestic vs Feral  1 0.3816603 2.336579 0.02530502    0.11       0.11    

beta_q1f$S %>%
  vegan::metaMDS(., trymax = 500, k = 2, trace=0) %>%
  vegan::scores() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(sample_metadata, by = join_by(sample == sample)) %>%
  group_by(origin) %>%
  mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
  mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
  ungroup() %>%
  ggplot(aes(x = NMDS1, y = NMDS2, color = origin, fill = origin)) +
    geom_point(size = 4) +
    scale_color_manual(values = origin_colors) +
    geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2), alpha = 0.9, show.legend = FALSE) +
    theme_classic() +
    theme(
      axis.text.x = element_text(size = 12),
      axis.text.y = element_text(size = 12),
      axis.title = element_text(size = 20, face = "bold"),
      axis.text = element_text(face = "bold", size = 18),
      panel.background = element_blank(),
      axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
      legend.text = element_text(size = 16),
      legend.title = element_text(size = 18),
      legend.position = "right", legend.box = "vertical"
    )